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Abstract - We study magnetic order induced by non-magnetic impurities in quantum paramag- 
nets with incommensurate host spin correlations. In contrast to the well-studied commensurate 
case where the defect-induced magnetism is spatially disordered but non- frustrated, the present 
problem combines strong disorder with frustration and, consequently, leads to spin-glass order. 
We discuss the crossover from strong randomness in the dilute limit to more conventional glass be- 
havior at larger doping, and numerically characterize the robust short-range order inherent to the 
spin-glass phase. We relate our findings to magnetic order in both BiCu^POe and YBa2Cu3C>6.6 
induced by Zn substitution. 
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Introduction. — Magnetic order induced by non- 
magnetic impurities is a remarkable phenomenon which 
has been observed in a wide variety of paramagnetic Mott 
insulators. It can be rationalized assuming dimer pair- 
ing of spins 1/2 in the host magnet: each non-magnetic 
impurity introduces one unpaired spin by breaking up a 
dimer. These impurity moments mutually interact by the 
exchange of gapped bulk excitations, with the interaction 
falling off exponentially as function of distance. If the host 
system features a bipartite lattice with antiferromagnetic 
spin correlations, the interactions between the impurity 
moments are non-frustrated, and consequently the ground 
state can be expected to be a long-range ordered antifer- 
romagnet. This has been verified both numerically [1-3] 
and experimentally, with Zn-doped CuGeC>3 [4,5] and Mg- 
doped TICUCI3 [6] being prime examples. 

Much less is known about impurity-induced magnetism 
in hosts with incommensurate (IC) correlations, realized 
e.g. in the Zn-doped spin ladder BiCu2POg [7-9] and also 
in the Zn-doped superconductor YBa2Cu306.6i where re- 
cent experiments [10] have demonstrated the appearance 
of static order at an IC wavevector close to that of the host 
correlations. Based on the available experimental data, it 
has been suggested that the physics is not very different 
from the commensurate case. 

In this paper, we instead argue that, although the mech- 
anism of moment formation by impurities is similar in 
the commensurate and IC cases, the low-energy physics 
of these moments is vastly different. The IC case leads 
to frustration of the inter-moment interaction and thus 



to spin-glass behavior. The actual ordering temperature 
where spins freeze into a glassy state is suppressed as com- 
pared to the commensurate case due to frustration, leading 
to a broad temperature regime of slow fluctuations. The 
glassy regime is characterized by pronounced short-range 
order at (or near) the wavevector imprinted by the host 
correlations, with a correlation length which vanishes as 
the impurity concentration x is tuned towards zero. Di- 
mensionality plays an important role: Glassiness is most 
pronounced if the host correlations are IC in three spatial 
directions, whereas a host with only one-dimensional (Id) 
incommensurability leads to a weakly glassy state with 
large magnetic correlation length. 

In the body of this paper, we present general arguments 
and extensive numerical simulation data, comparing the 
commensurate and IC cases, which lead to the above con- 
clusions. Our analysis sheds light onto a class of states 
which we believe is quite common to moderately disor- 
dered magnets, namely states with robust static short- 
range order combined with some amount of glassiness. We 
connect our results to the concrete experiments on doped 
spin-gap magnets in Refs. [7,10]. 

Effective model. — We shall assume that a mecha- 
nism of dimer breaking underlies the formation of mag- 
netic moments upon doping non-magnetic impurities. In 
its most general form, this mechanism relies on confine- 
ment of spinons in the host paramagnet and as such should 
apply to all spin-gap magnets with conventional spin-1 
(triplon) excitations [11]. 
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In the dilute limit of small x, there is a separation of 
energy scales, such that the impurity-induced moments 
provide the only degrees of freedom below the spin-gap 
energy A. Hence, the low-energy physics can be captured 
by an effective spin-S* Heisenberg model involving the im- 
purity moments Si at random positions 7"j only, 

W-efi = ~ JjjSj ■ Sj, (1) 
ij 

where S usually equals the bulk spin size. The inter- 
action Jij is dictated by the bulk magnetic properties. 
In a static linear-response approximation Jy = J 2 x{^ij) 
where Jo is of order of the bulk exchange coupling, and 
x(r) is Fourier transform of the static bulk susceptibility 
x(q)- Hence, \ Jij\ oc exp(— r^/^b) at long distances, where 
£b oc c/A is the bulk correlation length and c a mode ve- 
locity [12-15]. The separation of energy scales requires 
| Jij | <C A which implies I > £b where the average impu- 
rity distance £ = x~ x l d with d = 3. In this limit, the dis- 
tribution of energy scales in the model (1) becomes broad: 
the maximum coupling for each spin, Jj max = maxj |Jjj|, 
displays a strongly non-Gaussian distribution with large 
weight at small energies J" lax <C Jo [16]. 

Low-temperature order. — If the underlying lattice 
is bipartite and x{o) is peaked at the antiferromagnetic 
wavevector Q = (7r,7r, ...), then x(t) changes sign from 
site to site. All interactions can be satisfied classically by 
antiferromagnetic order of the Si , which has been shown 
to survive even in for S = 1/2 [2,3,7,17,18]. 

In contrast, if x(q) is peaked at an IC wavevector, the 
sign changes of x(^) are irregular, and the model (1) is 
generically frustrated [19]. Thus, spin-glass behavior is 
expected, which we will characterize below. 

It is worth noting that Eq. (1) is different from the 
Edwards- Anderson (E-A) spin-glass model [20]: The lat- 
ter describes spins on a regular lattice with random cou- 
plings, whereas Eq. (1) features spins at random positions 
with deterministic couplings [21]. As a result, the Jij in 
(1) display correlated disorder, in contrast to uncorrelated 
disorder usually assumed in the E-A model. Also, the 
broad distribution of couplings present in Eq. (1) is not a 
property of the standard E-A model. 

Monte-Carlo simulations. — Being interested in 
qualitative properties at low temperatures T including the 
ordered phase, we consider three-dimensional (3d) systems 
(with possibly anisotropic couplings) in the classical limit 
first, and turn to quantum effects later. We randomly 
distribute N[ mp unit vectors Si on a cubic lattice of lin- 
ear size L with periodic boundary conditions, yielding an 
impurity concentration x = Ni mp /L 3 . The effective in- 
teraction J^ is generated from Fourier transforming x(o) 
which is approximated by the inverse of the gapped bulk 
triplon dispersion e (q) [12]: 

J ^ J «V?7f (2) 



We perform classical Monte-Carlo (MC) simulations of 
the effective model (1), with single-site updates using a 
combination of the heat-bath and microcanonical (or over- 
relaxation) methods [22]. We consider typically 10 5 MC 
steps per spin for the measurements, after discarding 10 5 
MC steps for equilibration. For IC bulk correlations, 
we expect a (hardly-relaxing) spin-glass behavior, and to 
efficiently sample all spin configurations we employ the 
parallel-tempering algorithm [23]. Disorder averages are 
taken over N r \ impurity configurations, with A r i = 400 for 
L = 16 to N T i = 50 for L = 40. 

From the MC data, we calculate thermodynamic ob- 
servables, e.g., the specific heat C, as well as the neutron- 
scattering static structure factor 

iv imp ■ - x ' 

where (■ • • ) denotes MC average and [•••].„ average over 
disorder, and the sum runs over the A^ mp impurity sites. 
In a long-range ordered phase S(q)/Ni mp — )• M 2 5~q, 

where 8 is the Kronecker delta and Q is the ordering 
wavevector. Note that M — > 1 as T — > 0, because the clas- 
sical ground state is fluctuationless. The behavior of S(q) 
near Q allows to extract the correlation length £q charac- 
terizing the static impurity order (£q is distinct from £b 
[24]). To study signatures of spin-glass freezing, we also 
consider the spin-glass order parameter 

q a ^(q) = J-J2 S ? (1)s ? i2)eim > ( 4 ) 

where a and (3 are spin components, and W and ( 2 > denote 
two identical copies of the system ("replicas") containing 
the same impurity configuration. The spin-glass suscepti- 
bility is then defined as 

XSG(q) = N ilap ^2[(\q^(q)\ 2 )] av . (5) 

From xsg(<7) near q = we extract a spin-glass correlation 
length £sg- As q a '° acquires long-range order in a spin- 
glass state, a divergence of £sg signals the onset of spin- 
glass order, i.e., freezing. 

In our finite-size simulations, the ordering (or freez- 
ing) temperature T g is most efficiently extracted from 
the crossing points of £(T)/L data for different L (with 
£ = £q or £ = £sg)j according to the scaling law £/L = 
/(L 1 / iy (T — T g )), where f(x) is a scaling function and v is 
the correlation length exponent [25]. 

Spin-glass order with short-range correlations. 

— A detailed comparison between the cases of commen- 
surate and IC host correlations is in Fig. 1. We have 
used the bulk triplon dispersions e (q) in Figs, la and 
b, respectively, to generate the Jij according to Eq. (2). 
Fig. la has a single minimum at (7r,7r,7r) whereas Fig. lb 
has eight minima at (ir ± 6, n ± 6, n ± S) with S = -k/A 
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Figure 1: Classical MC results for the ordering process of im- 
purity moments described by H c g (1), in the cases of commen- 
surate (left) and incommensurate (right) bulk correlations, for 
x = 2% and different system sizes: L — 16 (O), 24 (□), 32 (A), 
and 40(0). ( a > b) Specific heat C as a function of T. In- 
sets: Host dispersion e (q), normalized by its bandwidth, along 
q = (q x ,iv,n) and q = (q x , §7r, §7r), respectively, (c, d) Mag- 
netic correlation length £q divided by L as function of T for 
Q — (-k.-k.tt) and (§7r, §7r, §tt), respectively, (e, f) T depen- 
dence of the spin-glass correlation length (,sg/L. The crossing 
point for different L defines the ordering/freezing temperature 
T s (vertical dashed lines). Insets in c, e, f show a zoom on £/L 
near T g . Error bars are smaller than the symbol size. 



[19]. In both cases we have used x = 2% and chosen 
e (q) such that £b — 3a. Energies are given in units of 

Jmax = \Jij(rij = 1)| [26]. 

All features in Figs. la,c,e consistently point towards 
the onset of long-range magnetic order at Tn « 0.82 J max : 
a unique crossing point in both £,q/L and £sg/L a t the 
same Tn [27] and a specific-heat maximum which both 
sharpens and moves towards Tn upon increasing L. 

The results for the IC case, Figs. lb,d,f, are strikingly 
different. A crossing is present only in ^sq/L, indicating 
a transition into a spin-glass state at a low T g « 0.07J max 
[21], while no crossing is observed in £,q/L. The peak in 
C(T) is broad and occurs at a temperature considerably 
larger than the freezing temperature (here T pea k ~ 2T g ). 
The origin is the emergence of short-range magnetic order 
far above T g [20], which is nicely visible in the gradual 
increase of £q, Fig. Id. The momentum dependence of 
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Figure 2: Structure factor S(q) for IC bulk correlations and the 
same parameters as in Fig. 1. a) Momentum dependence along 
q = §tt, §tt) at T = T a l% b) S(Q) at Q 
divided by iVi mp as function of T. 



the structure factor S(q), Fig. 2a, indeed shows peaks at 
the wavevectors Q corresponding to the minima of s (q). 
However, those peaks grow slower than the system size, 
Fig. 2b, indicating static short-range order with a vanish- 
ing magnetic order parameter M [28]. 

The fact that S(q) peaks at (or near) Q in the IC case 
is non-trivial given that the system is both strongly dis- 
ordered and frustrated. It may be rationalized by consid- 
ering the dense limit x — > 1: Here, the classical T-L e s (1) 
is minimized by a spiral state with ordering wavevector 
Q for which the Fourier transform of Jy has a maximum 
[29]. Our results show that, upon dilution, this spiral or- 
der then becomes both short-ranged and glassy. 

Doping dependence. — The model (1) has two char- 
acteristic length scales as input, namely t and £b, and 
one expects that their ratio £/£b is dominant in determin- 
ing the physical behavior. Clearly, £ 3> £b represents a 
strongly disordered regime, with a broad distribution of 
coupling constants and suppressed ordering temperature, 
whereas £ -c £b can be expected to lead to a more con- 
ventional state (still being glassy in the IC case). Many 
experiments are in an intermediate regime of £ > £b [7] . 

We have calculated the ordering/freezing temperature 
T g as function of doping x at fixed £b, with results shown 




Figure 3: a) Ordering temperature (rescaled by 0.25, □) for 
commensurate and freezing temperature (O) for IC host corre- 
lations, as function of the impurity concentration x. b) Mag- 
netic correlation length £q for IC host correlations, extrapo- 
lated to T = and L — oo, as function of x. (Dashed lines are 
guides to the eye, see text.) 
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in Fig. 3a. As expected, T s / J ma x is significantly smaller 
in the IC case as compared to the commensurate one 
due to frustration. T g (x) is approximately linear in this 
range of x for the commensurate case, while it shows sub- 
linear behavior for the IC case, which is well fitted by 
T g cx x 16 (dashed). For smaller £h also the commensurate 
case shows sublinear behavior (not shown). Experimen- 
tally, a linear T g (x) is frequently observed [7], while for 
very small x sublinear behavior has been found [5], con- 
sistent with our numerical results and previous ones for 
the commensurate case [7, 30] . 

Thus, T g (x) appears to be generically sublinear at small 
x before it crosses over to approximately linear behavior. 
It has been suggested that the typical value of the cou- 
plings Jij governs T g in the dilute limit x — > 0, while the 
average coupling is relevant in the more dense case where 
T g (x) is found to be linear [7,31]. We note, however, that 
at present there is no analytical understanding of the be- 
havior of T g (x) in the dilute limit. 

Fig. 3b shows the magnetic correlation length £q in the 
IC case, extrapolated to T = and L = oo. Remarkably, 
£,q{x) appears to vanish as x — >• 0, i.e., the spatial cor- 
relations between the impurity moments diminish in the 
dilute limit. Qualitatively, such a £q <C £ corresponds 
weakly correlated isolated impurities, whereas £q > £ can 
be associated with a spin glass formed from clusters of 
correlated impurities with cluster size £q. (For the pa- 
rameters in Fig. 3b, £q « I at x w 1%.) 

Dimensionality. — We have also considered spatially 
anisotropic cases where the IC host correlations only ex- 
tend along one or two directions (dubbed ldIC and 2dIC), 
whereas the correlations are commensurate (and weaker) 
in the remaining direction(s) - this situation applies to 
various quasi-ld or quasi-2d materials. 

The sample results in Fig. 4 demonstrate a clear trend 
with dimensionality of the IC correlations: The ldIC case 
in Fig. 4 (left) displays frustration in only one direction, 
therefore magnetic correlations are longer-ranged and con- 
sequently glassiness is weak. As can be seen from Fig. 4a, 
£q exceeds our largest system size at low T . However, 
from both the shift of the crossing points with system size 
(Fig. 4a), i.e., the absence of a unique crossing point, and 
the fact that S(Q)/Ni mp extrapolates to zero as L —> oo 
we can safely conclude that there is no conventional mag- 
netic order, but only spin-glass order at low T. Turning to 
the 2dIC case in Fig. 4 (right), we observe that this qual- 
itatively resembles the situation in 3d, where £q remains 
small below T g . 

We have studied the systematics of T g {x) in the spa- 
tially anisotropic cases, and find the general trend - sub- 
linear behavior at small x and linear behavior at larger x 
- to be obeyed here as well. Interestingly, T g (x)/J max is 
only suppressed by factors 1.3 — 4 by ldIC correlations as 
compared to the fully commensurate case (for correlation 
lengths as in Fig. 4a and 0.5% < x < 2%); this can be 
contrasted to the suppression by a factor of 10 or more 
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Figure 4: Ordering process of impurity moments, now for spa- 
tially anisotropic incommensurability, comparing ldIC with 
Q = (§7r,7r,7r) (left) and 2dIC with Q = (|tt, §tt, tt) (right). 
The bulk correlation lengths are £b/a ~ (4, 1.5, 1.5) (left) and 
£b/a ~ (3.3, 3.3, 1.3) (right) in the three directions. (a,b) Mag- 
netic correlation length f q (T) / L along an incommensurate di- 
rection for x — 0.5% and system sizes L — 24(D), 32(A), 
40(0), and 48(*). In (a), the arrows indicate the crossing 
points of the curves for different pairs of L. (c,d) Structure fac- 
tor S(Q)/Nim P . (e,f) Spin-glass correlation length £sg(T)/L. 
For both ldIC and 2dIC, a clear crossing point in £,sa/L is 
observed which allows to extract the T g (vertical dashed). For 
a discussion see text. 

found for the 3d incommensurate situation, Fig. 3a. 

Quantum effects and real-space renormalization 
group. — Our MC simulations presented so far neglect 
quantum fluctuations of the impurity moments. (Quan- 
tum effects in the bulk are crucial for the moment for- 
mation; those are taken into account in writing down 
Eq. (1).) A central question is whether the classical or- 
der of the Si survives for S = 1/2. For the commensurate 
case, both real-space renormalization group (RSRG) stud- 
ies of Eq. (1) [17,18] and numerical simulations of vacancy- 
doped quantum paramagnets [2,3,7] show that the answer 
is yes. In particular, the RSRG shows the generation of 
large spins in the course of the RG, which indicates mag- 
netic order at low T (instead of a random-singlet ground 
state) - the same has been found for random Heisenberg 
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models [32]. 

We have implemented the RSRG for the model in 
Eq. (1) with quantum spins 1/2 following Refs. [17,18,32]. 
In this procedure, we consider all possible pairs of spins in 
the system and calculate the energy gap between their 
ground state and the first excited multiplet. We then 
select the pair (Si,Sz) with the largest energy gap A12 
(which may then be identified with the system tempera- 
ture). One RG step is characterized by the decimation of 
the pair (Si, S2) according to the rules in Refs. [17,18,32], 
where the pair is either replaced by a single effective spin 
(representing a spin cluster) or removed completely (if the 
ground state of the pair was a singlet). We iterate the 
decimation procedure up to the last effective spin, or until 
we decimate the last spin singlet. 

Selected RSRG results for the model (1), Fig. 5, demon- 
strate the formation of large spins upon decimation for 
both commensurate and IC bulk correlations, very similar 
to the results in Refs. [17,18,32]. Both the average spin 
size S e g and the average number of original spins N e g of 
the remaining spin clusters diverge in the course of the 
RG, indicating magnetic order. This divergence can be 
described as S c g oc [32], and we estimate £ as 0.51 
for commensurate and 0.58 for IC bulk correlations, re- 
spectively. Importantly, the formation of a random-singlet 
state would instead imply that both S e ff and N e g remain 
small. The main difference between the commensurate 
and IC cases in Fig. 5 is that the high-temperature region 
where such random-singlet formation dominates is much 
extended for IC correlations, reflecting a broad regime of 
slow fluctuations above T g , as anticipated above. Finally 
we note that the RSRG procedure is not able to distin- 
guish between conventional long-range magnetic order and 
spin-glass order, because the effect of frustration - which 
requires at least three spins - cannot be fully captured 
within the pairwise decimation procedure. 

We conclude that the RSRG confirms the tendency of 
Eq. (1) toward low-temperature order: Quantum effects 
do not qualitatively change the predictions of the classical 
MC simulations, although they are important in a regime 
above T g where the formation of random singlets domi- 
nates in the dilute limit £ ;§> £t>- Hence, the classical MC 
simulations cannot be expected to describe the thermo- 
dynamics at and above T g , but they capture the physics 
of the low-T ordered state, with our T g being an upper 
bound to the true freezing temperature. 

Experiments. — BiCu2P06 is a spin-ladder material 
with ldIC correlations [8[. Hence, we expect weakly glassy 
behavior upon Zn doping [7, 9] . Indeed, deviations be- 
tween experimental data and the results of numerical sim- 
ulations (which do not account for frustration) have been 
observed [9] and tentatively assigned to the effect of mag- 
netic frustration. The fact that the experimental T g of 
BiCu 2 P06 is not very different from that of doped com- 
mensurate spin-gap magnets with similar energy scales 
and correlation lengths [7] is consistent with our finding 




A/A A/A 

max max 



Figure 5: RSRG results for the model (1) with quantum spins 
1/2, both for x = 1%, L = 50 (0,») and x = 2%, L = 40 
(□,■), with N T \ = 10 4 . Open (full) symbols correspond to 
commensurate (IC) bulk correlations, with parameters as in 
Fig. 1. (a) Size of the effective cluster spin moment S c s, av- 
eraged over the remaining (active) cluster spins, as a function 
of the energy scale A (equivalent to temperature), normalized 
by its maximal value, (b) Number of original spins- 1/2 per 
cluster, averaged over the remaining (active) cluster spins, as a 
function of A/A max . The dashed lines are power-law fits, and 
the plateaus at small A arise from the finite value of A?j mp . 

that ldIC correlations only cause a moderate suppression 
of T g . From our results we predict a strongly sublinear x 
dependence of T g for x < 1%. 

The superconductor YBa2Cu306.e displays a sizeable 
spin gap, with 2% of Zn substitution inducing static short- 
range magnetic order peaked at a 2dIC wavevector close 
to that of the host correlations [10]. For cuprates, experi- 
ments indicate that both IC magnetic correlations and IC 
magnetic order (the latter occurs at lower oxygen doping 
in Zn-free YBa2Cu306+5 as well as in other cuprates) can 
be described as collective magnetism (as opposed to sim- 
ple Fermi-surface nesting) originating from the tendency 
towards stripe formation [33]. Therefore it is reasonable 
to assume that the physics of impurity moments in spin- 
gapped cuprates can be captured by % ff (1) [11], and 
consequently we predict the magnetically ordered state in 
Zn-substituted YBa 2 Cu306.6 to be glassy. (In fact, Ni 
substitution in Lai.8sSro.i5Cu04 has been shown to in- 
duce spin-glass order even in the superconducting state 
[34].) Based on Fig. 3b we also predict the width of the 
magnetic Bragg peaks to increase for smaller x which can 
be checked by future neutron scattering experiments. 

Conclusions. — We have investigated the magnetism 
of impurity moments in incommensurate spin-gap mag- 
nets, which we have shown to be frustrated, leading to 
static order at low T akin to that of a cluster spin glass 
- in strong contrast to the commensurate case where the 
low-T state is a disordered, but unfrustrated antiferromag- 
net. 

Our results call for further experiments on impurity- 
doped magnets, in particular a characterization of spin- 
glass behavior of the low-T order, e.g., by ac suscepti- 
bility measurements. Theoretically, it will be interesting 
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to study the interplay of spin-glass and quantum effects 
beyond the real-space RG procedure. Particularly fasci- 
nating will be the connection of the zero-field spin glass to 
the Bose-glass phase expected in an applied held [35]. 

We thank P. Henelius, B. Keimer, C. Ruegg, and C. 
Timm for illuminating discussions. This research was sup- 
ported by the DFG through FOR 960 and GRK 1621. 
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